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ABSTRACT 

The multiple-planet systems discovered by the Kepler mission exhibit the following feature: planet 
pairs near first-order mean-motion resonances prefer orbits just outside the nominal resonance, while 
avoiding those just inside the resonance. We explore an extremely simple dynamical model for planet 
formation, in which planets grow in mass at a prescribed rate without orbital migration or dissipation. 
We develop an analytic version of this model for two-planet systems in two limiting cases: the planet 
mass grows quickly or slowly relative to the characteristic resonant libration time. In both cases the 
distribution of systems in period ratio develops a characteristic asymmetric peak-trough structure 
around the resonance, qualitatively similar to that observed in the Kepler sample. We verify this 
result with numerical integrations of the restricted three-body problem. We show that for the 3 : 2 
resonance, where the observed peak-trough structure is strongest, our simple model is consistent 
with the observations for a range of mean planet masses 2O-1OOA/0. This mass range is higher than 
expected, by at least a factor of three, from the few Kepler planets with measured masses, but part 
of this discrepancy could be due to oversimplifications in the dynamical model or uncertainties in the 
planetary mass-radius relation. 

Subject headings: planetary systems - planets and satellites: dynamical evolution and stability - 
planets and satellites: formation 



1. INTRODUCTION 

As of fall 2012, almost 100 multi-planet systems have 
been detected by ground-based radial-velocity observa- 
tions, and 365 transiting multi-planet systems have been 
detected by the Kepler spacecraft^. These large sam- 
ples enable statistical studies of correlations between the 
properties of members of multi-planet systems. One re- 
sult from these studies is that many planets appear to 
be close to mean-motion resonances, that is, their or- 
bital periods are close to the ratio of two small integers 
(see below for a more precise definition). In particular: 

• iLissauer et al.l ()2011bf ) comment that about 
one-third of the multi-planet systems studied 
by radial-velocity measurements contain near- 
resonant planet pairs, with about half of these near 
the 2 : 1 resonance. This is probably an underesti- 
mate because of a detection bias: there is an ap- 
proximate degeneracy between the signal from an 
interior planet near the 2 : 1 resonance and harmon- 
ics arising from non-zero eccentricity in the outer 
planet orbit, which makes it difficult to detect low- 
mass pla nets in resonance with exter ior compan- 
ions (e.g.. lAnglada-Escude et al.ll2010f ). 

• iFabrvckv et al.l ()2012l ) find that Kepler planet pairs 
with orbital period ratios within a few percent of 
2 : 1 or 3 : 2 are preferentially found just wide of the 
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resonance (i.e., period ratio slightly larger than 2 
or 1.5) and tend to avoid spacings just narrow of 
the resonance. 

• Near-resonances involving three or more planets 
are also present in the Kepler data. ILissauer et al.l 
(|2011bf ) find that the four planets in the system 
KOI-730 are in a chain of resonances, with pe- 
riod ratios 8:6:4:3 to within 0.1%. Lissauer et 
al. also find strong evidence for other resonances, 
including two three-planet resonances in t h e five - 
planet system KOI-500. IFabrvckv et~aT1 ()2012l ) 
find that three planets in the systems KOl-720 
and KOI-2086 have mean motions n that satisfy 
5ri2 — 3n4 — 2ni ~ and ni — In-i + ~ respec- 
tively, to within 0.01% (the planets are labeled in 
order of increasing period) . 

• In the GJ 876 system, planets 2, 3, and 4 have 
periods close to the ratio 1:2:4; planets 2 and 3 
are close to a secular resonanc e; and there a re other 
probable near-resonances (see lBaluevllMTI and ref- 
erences therein). Because of these near-resonances, 
gravitational interactions between the planets are 
detectable in the radial-velocity data, and these 
can be used to constrain the mutual inclinations 
and determine whether various critical arguments 
librate or circulate. 

• Within the solar system, Jupiter and Saturn are 
within 1% of a 5:2 resonance (the "great inequal- 
ity"), Uranus and Neptune are within 2% of a 2 : 1 
resonance, and Pluto is in a 3 : 2 resonance with 
Neptune. There are also many resonances among 
the satellites of Jupiter and Saturn. Another no- 
table near-resonance occurs between the two outer 
planets of the pulsar PSR B1257+12, which are 
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within 2% of a 3 : 2 resonance. This near-resonance 
produces gravitational interactions large enough to 
be easily detectable, which allowed the existence 
of the planets to be confirmed and their inclina- 
tions to be measured shortly after their discovery 
(|Rasio et alJll992l : iMalhotra et"aII[T992T ). 

Loosely speaking, an orbital resonance between two 
planets occurs when their mean motions or orbital fre- 
quencies ni,ri2 are nearly commensurate, i.e., nxjni is 
close to a ratio of small integers, p : p+q where and 
q > 0. The case q = is sometimes called a corotation 
or co-orbital resonance; examples in the solar system in- 
clude the Trojan asteroids (1:1 resonance with Jupiter) 
and the Saturnian satellites Janus and Epimetheus, but 
no extrasolar co-orbital resonances are known. When 
q > it is called the order of the resonance, since for 
planets on nearly circular, coplanar orbits the strength 
of the resonance potential is proportional to e q or I q 
where e and I are the eccentricities and inclinations of 
the resonant planets. Inclination resonances occur only 
for even q. 

In a resonant configuration, the longitude of the plan- 
ets at every qth conjunction librates slowly about a di- 
rection determined by the lines of apsides and nodes of 
the planetary orbits. In the action-angle variables for the 
Keplerian potential, this geometry is naturally described 
by the libration of a so-called critical argum ent which 
is a l i near combination of the angle v ariables (jMalhotral 
[1991 IT991 [Murray fc Dermottill999l) . 

Definition of resonance — Generally, a planet pair is said 
to be "in resonance" if some dynamically significant crit- 
ical argument librates. If the critical argument circulates 
the planet pair is said to be near but not in resonance. 
This definition h as several shortcomings for our purposes 
(see f or example iHenrard fe Lemaitrelll983T : lDelisle et al.l 
2(1 1 2;-: (i) it is not always consistent with the expec- 
tation that resonant planets have orbital periods close 
to the ratio of two small integers, since planets can be 
resonant at any period raticQ; (ii) by suitable canonical 
transformations one can change the appropriate critical 
argument from libration to circulation, so this definition 
of resonance is coordinate-dependent; (hi) interaction of 
nearby resonances arising from the degenerate frequen- 
cies in the Kepler problem can cause a critical argument 
to jump between libration and circulation at irregular 
intervals. 

Given these comments, it is nugatory to discuss 
whether Kepler's multi-planet systems are "resonant" 
or "non-resonant" , and in this paper we shall avoid the 
term "resonant" in favor of the looser description "near- 
resonant" . 

Resonances and migration — Convergent migration — 
evolution of the semi-major axis of one or both planets 
such that the period ratio approaches unity — can lead 
to permanent capture into resonance, and is b elieved to 
be th e cause of the Neptune-Pluto resonance (jMalhotral 
119931 ) and the resonances between satellites of Jupiter 

5 As the eccentricity approaches zero, the apsidal precession rate 
due to a perturber grows without limit. Thus a critical angle whose 
time derivative involves the mean motions and precession rate can 
librate even when the mean motions are far from commensurability. 



and Saturn (although in these cases the migration is out- 
ward, whereas migration in most cxoplanct systems is 
believed to be inward). 

The existence of near-resonant plan et pairs is often 
ascrib ed to convergent migration (e.g., ISnellgrove et al.l 
l200l . but there are problems with this hypothesis: 

• The fr action of plane t s in resonance is quite 
small: iFabrvckv et al.l ()2012f ) find peaks at the 
3 : 2 and 2 : 1 resonances of ~ 20 planet pairs 
each, out of a total sample of ~ 750 planet 
pairs; this small fraction is noteworthy since cap- 
ture into the 2 : 1 resonance is certain during con- 
vergent migration if the planets cross the reso- 
nance slowly enough and their initial eccentrici- 
ties and inclinations are small enough. "Slowly 
enough" means a migration timescale larger than 
10 4 yr(P pl /100d)(10M e /TO pl ) 4 / 3 where m p] and 
P p i are the migrating planet's mass and orbital 
period; "small enough" means initial eccentricity 
smaller than O.O5(m p i/lOM0) 1 / 3 (see Appendix). 
Typical Type I migration times for Kepler planets 
i n a low-mass prot oplanetary disk are 10 3 -10 4 yr 
(jTanaka et al.ll2002f ). in which case migrating plan- 
ets might avoid resonance capture according t o this 
criterion, but orbit integrations bv lReinl (|2012|) that 
include eccentricity damping show that most Ke- 
pler planets are captured even if the migration time 
is as short as 10 3 yr. 

• Convergent migration typically leads to capture 
into a 2 : 1 resonance. Capture into resonances 
with smaller separations, such as 3:2 or 4:3, be- 
comes increasingly difficult: for example, to cap- 
ture a planet into the 3 : 2 resonance requires cither 
that the initial conditions are fine-tuned so that the 
planets form in the narrow interval between the 2 : 1 
and 3 : 2 resonance (only 13% of the outer planet's 
semi-major axis), or that the migration rate is fast 
enough that the planet jumps the w eaker 2 : 1 and is 
captu red at the 3 : 2 resonance. Yet IFabrvckv et all 
(|2012[ ) find that the excess of planet pairs near the 
3 : 2 resonance is at leas t as strong as the excess at 
the 2 : 1 resonance (see IRein eFatl 12011 for a de- 
tailed numerical study of these issues). 

• In the s ample of Kepler multi- planet systems exam- 
ined by IFabrvckv et all (|2012[ ). the excess of planet 
pairs at period ratios just larger than 2 is accompa- 
nied by a deficit at period ratios just smaller than 
2. Both the peak and the trough have an equivalent 
width of about 20 planets. This strongly suggests 
that the features at this resonance arise from rear- 
ranging the periods of planets near the resonance, 
rather than by capturing planets at the resonance, 
which should produce a peak but no trough. 

IBatygin fc Morbidellil (120121 ) and iLithwick fc Wul 
(|2012f ) have pointed out that dissipation due to tides 
from the host star or the protoplanetary disk tends to 
repel near-resonant planets, in the sense t hat their pe- 
riod r atios evolve away from unity (see also iDelisle et al.l 
I20I2T) . Thus, if migration is common, dissipation could 
explain why there are so few near-resonant planet pairs; 
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and if there is little or no migration, dissipation could 
explain why planet pairs that initially happen to lie near 
a resonance are now preferentially found wide of the 
resonance. This hypothesis is discussed further in §7.31 

Although these and other analyses in the literature 
shed considerable light on the behavior of near-resonant 
planets in the presence of dissipation and/or migration, 
our view is that they put the cart before the horse: the 
question that should be addressed hrst is, what is the dis- 
tribution of mean motions or period ratios expected near 
a resonance in the absence of dissipation or migration? 

Efforts to address this question have a long and rich 
history, mostly in the context of the asteroid belt, in 
which the distribution of mean motions contains gaps at 
the 4:1,3:1, 5:2, 7:3, and 2 : 1 resonances with Jupiter 
(the Kirkwood gaps) and peaks at the 3 : 2 and 1 : 1 res- 
onances (the Hilda family and the Trojan asteroids). 
Most of these features can be largely explained through 
chaotic evolution of the asteroid orbits on timescales 
as long as 1 Gyr, lon g after the formation of the solar 
system was complete (lWisdomlll983t iMurrav fc Holmanl 
119971 : iLecar et al.1 120011 ). Although the effects of dissi- 
pation and migration are discernible in the orbital dis- 
tributi on of asteroids, these features are comparatively 
subtl e (|Liou fc Malhotrall997l : lMinton fc Malhotral20"09l 
l2010f ). Thus the asteroid belt demonstrates that dissipa- 
tion and large-scale migration are not essential to pro- 
duce near-resonant features in the mean-motion distri- 
bution. Unfortunately for our purposes, the chaotic evo- 
lution of the asteroids appears to depend strongly on the 
details of the planetary configuration in the solar system 
and thus the insights gained from studies of the asteroid 
belt cannot be immediately applied to exoplanets. 

In this paper we focus on a simple and highly idealized 
model for the formation of near-resonant features: at an 
initial time t = we place a single planet and a large 
number of test particles on circular orbits around the 
host star. The test particles are interior to the planet and 
smoothly distributed in semi-major axis (see Equation 
I41j) . The mass of the planet is m p \h(t) where h(t) ramps 
up from zero at time zero to unity at large times; we 
focus on the two limiting cases in which the planet mass 
grows rapidly — h(t) is a step function at t = — and 
slowly compared to the characteristic orbital and secular 
frequencies of the test particles; see Equation (|38l) . 

2. ANALYTIC RESULTS 
2.1. Orbital elements 

We follow the motion of the test particles using oscu- 
lating Keplerian orbital elements, with a,e,u;,£ being 
the semi-major axis, eccentricity, longitude of periapsis 
and mean anomaly. For canonical elements we use the 
modified Delaunay variables A = I + lo,j = —a;, A = 
M 1/2 ,r = (/ia) 1/2 [l - (1 - e 2 ) 1 / 2 ], where fi = Gm* is 
the gravitational mass of the star. 

The Hamiltonian for the restricted three-body problem 
(star, planet, test particle) is 



H (A, 7, A, T, t) = - — - /i p i 



1 



r — r 



Pil 



r ■ r p i 

r pi 



where /i p i = Gm p \ is the gravitational mass of the planet 
and r, r p i are the position vectors of the test particle and 
planet relative to the star. The first term describes the 
unperturbed Keplerian motion of the test particle about 
the Sun, and the remainder describes the perturbation 
from the planet. 

Transit observations of exoplanets measure the time 
interval Ptr between successive transits, that is, the in- 
terval in which the longitude increases by 2tt. We must 
relate the transit-based mean motion n tr = 2-7r/Ptr to the 
osculating mean motion n, given by n 2 = /i/a 3 . Write 
the perturbing Hamiltonian as 

ff p i(A J 7,A,r,t) = ff lp ( 7 ,A ! r ) t) + F sp (A,7 J A,r j t) (2) 

where 

H\ p = (H p i)t 7 H sp = H p i — Hip, (3) 

with (•)( denoting a time average over an interval compa- 
rable to the total observational timespan; H sp and H\p 
are respectively the short-period and long-period parts 
of the Hamiltonian. Then n t r is the time average of the 
rate of change of the mean longitude, n tr = (d\/dt) t , and 
we have 



_/dX\ _/dH\ _l^_ dHip_ 

™ tr " \ dt ) t ~ \dK/ t ~v + ~dA~~ 



dH 



ip 



dA 



(4) 



where n is the osculating mean motion. 

For the resonances we shall examine, the second term 
in these equations is small; thus we may neglect the dis- 
tinction between n and n tr - In other words, the distri- 
bution of transit-based period ratios near resonances is 
determined mainly by the distribution of the osculating 
mean motions rather than by any differences between the 
osculating and transit-based mean motions. 

2.2. Resonant dynamics 

We consider the vicinity of a first-order mean-motion res- 
onance, where the ratio of the mean motions of the planet 
and the test particle is p : (p+1). Thus p > corresponds 
to an interior resonance (test particle inside the planet) 
and p < — 1 is an exterior resonance. 

We now assume that the resonances arc sufficiently well 
separated compared to the size of the perturbations in- 
duced by the planet that we can ignore all of the per- 
turbations except those associated with the p : (p+1) 
resonance (we call this the single-resonance approxima- 
tion). It is useful to make a canonical transformation to 
slow and fast variables, 



</>=(p + l)A p i -pA + 7, 
ip = \- A p i, 



$ = T; 

® = A + pT, (5) 

where A p i = n p \(t — to) is the mean longitude of the 
planet, and n p i is its mean motion. Then the new Hamil- 
tonian is 



H =npi[(p + 1)$ - tf] - 



= -^3+/* P i.Hpi(A ) 7,A ) r,t) (1) 



2(* -p$) 2 
+ ^ p i-Hpi(</>, a p i, epi) (6) 

where a p \ is the semi-major axis of the planet's orbit and 
Hp\ represents the planetary perturbation. For simplicity 
we have assumed that the planet is on a circular orbit 
and that the test-particle and planet orbits are copla- 
nar. Since ip is a fast variable, we will drop ^-dependent 



4 



Petrovich, Malhotra & Trcmainc 



terms, which is equivalent to replacing _ff p i by -ffi p ; con- 
sequently, the resonance Hamiltonian is independent of 
ijj and ^ is a constant of the motion. We shall use the 
notation 

a c = — ; (7) 
M 



EL 

$3' 



n c and a c arc constants of the motion which equal the 
osculating mean motion and semi-major axis of the test 
particle when its eccentricity is zero. 
If the test-particle orbit is nearly circular, $ ~ ^y/JIae 2 



with a few 

H ICS = [(p + l)n p i - pn c ]$ + 2 + eV2$ cos 0, (8) 
where we have dropped an inessential constant, and 

ip 2 n c ^ _ /ipi f p 
2# ' ° ~ a p i v^" 



is small, and wc can approximate Equation 
terms in an expansion in powers of 



(9) 



Note that \& ~ ^/Jia(l + jpe 2 ), and since the eccentricity 
is small, we have ^ > and $ < in all cases of interest . 
The coefficient f p is given by ([Murray fc Dermottlll999| ) 

/ p = -(p+l + |£))6^ 1) (a) ) a=[p/(p+l)] 2 / 3 (10) 

when p > (planet exterior to test particle) , and 



f P = -a(p+f)-j)D)b«J+ ll) (a) 



Op -2 

2a 



a = [{p+l)/p] 2 ' 3 
(11) 

when p < — 1 (planet interior to test particle). Here 5ij 
is the Kronccker delta function, D = d/dloga, and the 
Laplace coefficient 



(«) = 



1 



cos mx dx 



(1 - 2a cos a; + a 2 ) 1 / 2 



(12) 



"1/2 ■ • 7T 

We have 

/i = -1.1905, f 2 =- 2.0252, / 3 = 
/_ 2 = 0.26987, /_ 3 = 1.8957, /_ 4 

Following IHenrard fc L cmaitrc (|1983D . we define a di 
mcnsionlcss time and canonical momentum, 

1/3 



-2.8404; 
2.7103. 



(13) 



(3e 2 



t, 



R 



e 



2/3 



and a modified canonical coordinate r 

r = —<b 



■IT 



if e > 
if e <0. 



(14) 



(15) 



The new, dimensionless Hamiltonian in the canonical 
variables (r, R) is then given by 

(16) 



K = -3AR + R 2 - 2V2R cos r 
where the dimensionless resonance distance A is 
(j> + 1) 77.pi - pn c 



A 



|27/3e 2 /4| 1 /3 



(17) 



A few notes: 



• The strength of the resonance can be parametrized 
by the dimensionless ratio 



|27/3e 2 /4| 



1/3 



pn c 



sgn(p) 



9/pTOpldc 



2/3 



(18) 



Note that s is positive for interior resonances and 
negative for exterior ones. 

• For small eccentricity, we can write 2R ~ (e/s e ) , 
where the eccentricity scale is 



1 



f p m v \a 



e 
2P 

1/3 



1/3 



3p 2 m*a p i 



Note that 



9p 



-s. 



(19) 
(20) 



• As the eccentricity e — > and the planet mass 
m pi — > 0, n c = n is the unperturbed mean mo- 
tion, and the "exact resonance" condition pn c = 
(p + l)n p ] corresponds to A = 0. 

• The topology of the phase space determined by this 
dimensionless resonant Hamiltonian depends only 
upon the value of A. 

• The range of the resonant perturbation is |(p + 
1)%>1 — pn c \ ~ pn c s. This means that the "width" 



of the resonance is proportional to m 



2/3 

pi 



The phase-space trajectories follow level curves of the 
dimensionless resonant Hamiltonian K (Equation I16[) . 
Figure Q] shows plots of the level curves for various val- 
ues of A to illustrate the phase-space topology. In 
these plots, we use the Cartesian variables (x,y) = 
V2i?(cos r, sin r) , which are also canonical (x is the mo- 
mentum and y is the coordinate). Thus, the origin in 
these plots corresponds to zero eccentricity, and the dis- 
tance from the origin is e/s e . The phase-space structure 
is simple when |A| 3> 1: the trajectories are nearly cir- 
cles centered close to the origin. For A < 1, there is 
only one fixed point and no homoclinic trajectory, but 
for A > 1 there are three fixed points and a homoclinic 
trajectory exists. All the fixed points are on the a;-axis; 
they are given by the solutions of dK/ dx = which are 
the real roots of the cubic equation 

x 3 - 3Ax - 2 = 0. (21) 

Figure [2] plots the locations of the real roots as a function 
of A. We shall be using only the branches shown by solid 
lines: for A < 



an (A) = (1+ v 7 !- A3) 1 /3 + A (i +v /i_ A3)- 1 / 3 , (22) 
and for A > 

x 2 (A) = -2VAcos[((9 - 2tt)/3], 6 = cos^-A" 1 / 3 ). 

(23) 
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Fig. 1. — Level curves of the dimcnsionlcss resonance Hamiltonian, Equation 116> . for various values of the resonance distance A, 
Equation {17}. The coordinates are (x,y) = V2R(cosr, sinr). 



The level surface of the Hamiltonian that passes through 
the X3 homoclinic point is a classical quartic curve called 
the limagon of Pascal. 

3. EVOLUTION OF INITIALLY CIRCULAR ORBITS IF THE 
PLANET MASS GROWS QUICKLY 

In this section we determine the distribution of test- 
particle orbits that arises if the planet appears quickly. 
Here "quickly" means over a timescale long compared to 
the orbital time but short compared to the inverse of the 
frequency scale sn c (Equation I18j) : for reference, at the 
2:1 resonance of Jupiter (sn c ) _1 ~ 100 yr. We recog- 
nize that this is probably an unrealistic model of how 



most planets form, but we present it as a foil to the slow 
planetary growth described in the following section. 

Assume that the planet mass grows suddenly from to 
m p i at time r = 0. Since the test particles have zero ec- 
centricity at this instant, R = so the resonance Hamil- 
tonian K = and the fast action $ = A = (//at) 1 / 2 
where is the initial semi-major axis. Then n c is a 
constant of the motion that is equal to = (/i/af) 1 / 2 , 
the mean motion just before r = (Equation . The 
resonant mean motion is n rcs = (p + l)n p \/p and the 
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Fig . 2. — Fixed points of the dimensionless Hamiltonian, Equa- 
tion JTB}, as a function of the resonance distance A. The fixed 
points have r = 0, n and the variable plotted is x = v2~f?cosr. 
The arrows indicate the evolution of initially circular orbits dur- 
ing the slow growth of the planet, while the dashed lines represent 
fixed points that are not populated by test particles during this 
process. 

dimensionless resonance distance (Equation [TT]) is 



A = 



(24) 



Thus a uniform initial distribution in mean motion im- 
plies a uniform distribution in A. 

Observers more commonly work with the period ratio 
V, defined to be greater than unity, i.e., V = n/n p \ for 
interior resonances and n p i/n for exterior ones. Then 



A = 



1 /Pros 

S \ Vi 

1 (Vi_ 

s \V res 



- 1 



1 



Vrcs = (P+1)/P, P> 

V ICS =p/(p+l), P < -1. (25) 



For \s\ <§; 1, the effects of the resonance are only impor- 
tant if \T /V les — 1| <C 1 and in this case we can combine 
the two formulae above as 



A=ifl-^L 

\s\ V r rci 



(26) 



For t > the particles evolve at fixed energy K and 
action ^ . For initially circular orbits, = 0, ^ = ^/Jlal, 
and the constancy of W implies that a[l + p — p(l — 
e 2 ) 1 / 2 ] 2 =a . and 



rii [1 + p - p 



=2l3 



ni (l + lpe 2 ). (27) 



Thus, the final mean motion or semi-major axis is de- 
termined by the final eccentricity e/, and our problem is 
reduced to determining the distribution of eccentricities 
e excited by the planet. Using Equation (|20j) . Equation 



o 
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Fig. 3. — Distribution of period ratios in the vicinity of a res- 
onance when the planet mass grows quickly, as described in J3] 
Period ra tios are plotted in units of the resonance strength |s| 
(Equation 1181 . 



(|2"T|) can be simplified to 



(28) 



Now suppose a particle initially has mean motion sep- 
arated from the resonance by a fractional distance fs 
where / is of order unity, i.e., rii — n rcs = fsn.i. Then 
from Equation (pM)) , A = — / so the resonance Hamilto- 
nian K (eq. [16]) is completely specified. The initial con- 
ditions are x = V 2i? cosr = 0, y = \ / 2Rsinr = 0, and 
we integrate the equations of motion dx/dr = —dK/dy, 
dy/dr = dK/dx to the present time; then the present 
value of R determines the present mean motion in units 
of s through Equation ([2"8"]l. Thus the steady-state distri- 
bution of mean motions is the same for all resonances of 
this type, independent of the planet mass, the resonance 
integer p, whether the resonance is interior or exterior, 
etc., so long as the mean motions are expressed in terms 
of the dimensionless parameter s of Equation (|18|). 

Using this procedure we can compute the expected dis- 
tribution of mean motions near a resonance in the limit- 
ing case where the planet mass grows fast. This has been 
done in Figured where we have chosen 10 7 particles with 
zero initial eccentricity and initial mean motion chosen 
randomly in the interval IPi/Pres — 1| < 3|s|. Each parti- 
cle has been followed for a time chosen uniformly random 
between 50 and 150. 

We observe from Figure [3] that the distribution of par- 
ticles exhibits a peak at period ratios larger than reso- 
nance, and a trough at period ratios smaller than reso- 
nance (this statement holds for both interior and exte- 
rior resonances). The shape resembles a P-Cygni profile 
in the context of spectroscopic emission lines. Using this 
analogy we can characterize the redistribution of systems 
close to resonances by computing the equivalent widths 
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(EWs) in these regions as 



/ 


~yi(P) J 


IV > V rcB 


. Vi 



where yi is the initial number of particles per unit pe- 
riod ratio (assumed uniform) and yf{V) is the final dis- 
tribution. Here, the initial distribution defines our back- 
ground or "continuum" . In our convention a trough ( "ab- 
sorption line") has negative EW and a peak ("emission 
line") has positive EW. If particles are rearranged near 
the resonance but do not migrate or escape the sum of 
the EWs should be zero. The results from the simulation 
in Figure |3] give 

EW+ = -EW_ = 0.685|s|"P rcs . (30) 

4. EVOLUTION OF INITIALLY CIRCULAR ORBITS IF THE 
PLANET MASS GROWS SLOWLY 

We now examine how initially circular orbits of test parti- 
cles in the vicinity of resonance evolve as the planet mass 
grows gradually from small values. The slowly varying 

parameter is the resonance distance A oc m p ^ 3 . For 
very small values of planet mass, |A| is large [for all but 
the measure zero case in which (p+ l)^ P i = pn c \- As the 
planet mass grows to its final value, |A[ decreases. 

4.1. Analysis 

Provided that there is no scparatrix-crossing event, an 
initially circular orbit of a test particle will evolve with 
two adiabatic invariants: = ^Jm{\ + p — pyl — e 2 ) 
and A = § Rdr = § ydx, the area enclosed by the phase 
trajectory in the (x, y) plane. If a separatrix crossing 
occurs, there is a discontinuous change in A; this does 
occur for a range of parameters, as detailed below. 

In the discussion below, we use subscripts i and / to 
denote initial and final values of parameters. For initially 
circular orbits, the adiabatic invariant ^ = y/JuXi and 
the constant n c (Equation [7|) is equal to the initial mean 
motion rij. 

Consider first the implications of the adiabatic invari- 
ant ^P. The relations (|24| and (|28| continue to hold, ex- 
cept that now the resonance strength s is growing slowly 
with time. Thus, the final mean motion or semi-major 
axis is determined by the final eccentricity e / or its scaled 
version Rf, and our problem is reduced to determining 
Rf. We now show that it is possible to determine Rf 
from the adiabatic invariance of A for most (but not all) 
initially circular orbits. 

For initially circular orbits, Ai — 0. Therefore, these 
orbits will adiabatically follow the evolution of the fixed 
point close to the origin. There are two cases to consider, 
one for each "side" of the resonance: 

Test particle is initially wide of the resonance — This means 
that the test particle is further from the planet than the 
resonance, that is, for an interior resonance (p > 0) the 
test particle is inside the resonance and rii > (l + l/p)n p i, 
while for an exterior resonance (p < —1) the test particle 
is outside the resonance and rii < (1 + l/p)n p i. For 
both cases Aj is negative, and large in absolute value. 
Initially circular orbits are near the fixed point x\ in 
Figure [2] As A increases, these orbits evolve along the 
fixed point x\, i.e., upward along the curve labeled x\ 



in the fourth quadrant of Figure [2J as indicated by the 
upward pointing arrow. No matter how large the planet 
mass grows, A remains negative, so the test particle can 
never evolve onto the dotted portion of the x\ curve. 
(The point marked with a black circle at A = 0, x\ = 2 1 / 3 
corresponds to the measure-zero case in which rii = (1 + 
l/p)n p i.) Thus, the final orbits have a fixed eccentricity 
(Equation 

e/ = s e ar 1 (A / ) or R f = |x? (A/), (31) 

and a final mean motion given by Equation (|28D . 

Test particle is initially narrow of the resonance — In this 
case Ai is large and positive and initially circular orbits 
are near the fixed point cc 2 in Figure [5] As A decreases, 
these orbits evolve along the fixed point X2, i.e., down- 
ward along the inner curve in the second quadrant in 
Figure [2 as indicated by the downward pointing arrow. 

No matter how large the planet mass grows, A remains 
positive. This branch represents the final orbits for A/ > 
1 (equivalently, pn,; < (p + l)n p \ — psn c ), and for these 
orbits the final eccentricity is fixed at (Equation |2"3"|) 

e f = s e x 2 {A f ) or R f = |xf(A/). (32) 

A complication in this case is that the fixed point x 2 
vanishes at A = 1, i.e., for prii — (p + l)n P i — psn c 
where s is the resonance strength defined by (|18|) . For 
< Af < 1 [i.e., (p+l)n p i-psn c < pn t < (p+l)n p i], the 
adiabatic evolution leads the trajectory to coincide with 
the separatrix when A = 1. There ensues a discontinuous 
increase in the value of A, which jumps from A = Ai = 
to A = 6tt (which is the area enclosed by the separatrix 
at A = 1). After this jump A is again an adiabatic 
invariant so its final value is Af = 67r. The final orbits 
in these cases are not described by a stationary value 
of the eccentricity, because the trajectory in the (x, y) 
plane is neither a fixed point nor a circle centered at the 
origin. Similarly, the final mean motion is not fixed. An 
approximate estimate of the mean final eccentricity is 

(e)/ = s e JAf/-K = VGse, (33) 

but to determine the final distribution of eccentricities 
and mean motions in this region it is simpler to integrate 
the equations of motion in the resonance Hamiltonian K 
(see the following subsection). 

Combining Equations (US]), (|5i). and ([3"3)l. the final 
mean motion of particles that do not cross the separatrix 
is 

n =n,(l + \sx\[(n rcs ~ nij/am]), p(n rcs - n,) > 0, 

=m(l + ^sxl[(n rcs - m)/sni]), 

p[n Ies - (1 + s) ni \ > 0, (34) 

where n, is the initial mean motion, and the character- 
istic resonance strength s (Equation [T81 is evaluated for 
the final planet mass m p \. If the initial distribution in 
mean motion is uniform, the final density of the mean- 
motion distribution is given by 

dN(n) oc drii = ^^-dn, (35) 
dn 

which is easily evaluated numerically from Equation 
([34])— see Figure g] 
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(V/V res -l)/\s\ 

Fig. 4. — Distribution of period ratios in the vicinity of a reso- 
nance when the planet mass grows slowly, as described in Jl] Pe- 
riod ratios are plotted in units of the resonance strength |s| (Equa- 
tion [THJ. The red dashed curves show the expected distribution for 
part icles that do not cross the separatrix, as given by Equations 
(1341 and (l35l ■ 



For these particles there is a gap in the distribution of 
final mean motions, 



«-rcs(l 

n ies (l - 



3 

2 2/3 



|s|) < n < n rcs (l + 2-3— p,s>0, 
|s|) <n<n res (l + ||s|) p,s<0. (36) 



Note that the gap is not symmetric about the exact res- 
onance value, n rcs = (p + l)n p i/p. The size of the gap 
is 



An gap = 2.587 (p+1) 



m p i a c f P 



TO* Cipip^ 



n p \. 



(37) 



The particles that cross the separatrix have librating ec- 
centricities and mean motions; as we show below, these 
particles partially fill the gap and also broaden the peak 
wide of resonance. 

4.2. Numerics 

For our numerical experiments, we assume that the 
planet mass varies as 



to p i(£) = rripj tanhi/tpi. 



(38) 



The timescale t p \ can be thought of as the formation time 
of the planet. In simulations with t p \ = the planet 
forms suddenly, i.e., the integration is started with the 
final planet mass m p \. 

We follow the evolution induced by the resonance 
Hamiltonian K for a large number of test particles on 
initially circular orbits, uniformly distributed in mean 
motion. We start the integrations at a time when the 
planet mass to p i^ = 10 _ m p ij, which means that the 



any particle the final resonance distance A/ = 10" 4 Aj; 
we shall follow particles with |A/| < 3 which implies that 
the initial distribution should be chosen uniform in A.; 
between ±3 x 10 4 . We use a growth time t p \ = 100 and 
follow the particles for a time chosen uniformly random 
between 250 and 750; we have checked that the results 
are insensitive to these choices. The resulting distribu- 
tion of mean motions is shown in Figure [4] as expected, 
the distribution mostly agrees with the distribution de- 
rived analytically in the preceding subsection, but the 
particles that have crossed the separatrix partially fill in 
the resonance gap (Equation I36|) and enhance the peak 
to the right of it . 
The equivalent widths arc 



EWj 



-EW_ = 0.956|s|P re 



(39) 



almost 50% larger than the result when the planet grows 
fast (Equation [30]) . 

5. NUMERICAL INTEGRATIONS OF THE RESTRICTED 
THREE-BODY PROBLEM 

We ran numerical experiments to follow the motion of 
test particles subject to gravitational forces from the cen- 
tral star and an orbiting massive planet. We work in the 
astrocentric reference frame, in which the equations of 
motion for a test particle can be written as 



1 d 2 x 



x-x p i(t) x p i(t) 



\x-x pl (t)\ 3 |z P i(i)P 



(40) 



where fix = m*/(TO* + m p \) and /12 = 1 — [J-i = 
m P \/(m* +to p i). As usual a p \ and n^ x = G(TO*+m p i)/ap[ 
are the semi-major axis and squared orbital frequency of 
the planet. The positions of the test particle x and the 
massive planet x p \ are normalized by a p i. Times are given 
in planet years, 27r/n p i. 

In these simulations the planet mass m p i or /i2 is ini- 
tially zero and grows during the simulation according to 
the formula (|38|) . During this growth the planetary semi- 
major axis dpi is kept constant, although there are other 
plausible choices, e.g., a p \ cx (to* + m p i) _1 , as would be 
expected if the planet gained mass isotropically. 

If the planet eccentricity e p i is zero and the growth of 
the planet mass is sufficiently fast or slow, the results of 
these simulations should be directly comparable to the 
analytic results we obtained in an d SjU 

5.1. Initial conditions 
We typically consider planet masses m p i in the range 



io- 4 -io- 



or 0.1-1 Jupiter mass for a solar-mass 



initial characteristic frequency Si 



10" 



Thus for 



host star. For comparison most Kepler planets have 
masses between 0.01 and 0.1 Jupiter masses, and most 
planets discovered by radial-velocity measurements have 
masses between 0.1 and 10 Jupiter masses. 

We assume for simplicity that the test particle is in- 
side the massive planet, i.e., we consider interior reso- 
nances only. This is a plausible simplification: inner 
planets are usually smaller than outer planets in the 
Kepler and radial-velocity samples because planets with 
smaller semi-major axes are easier to detect. 

The semi-major axes a of the test particle and planet 
are determined using a fitting function for the prob- 
ability distribution of semi-major axes in the Kepler 
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sample, after accounting for geometric selection effects 
(jTremaine fc Dondl20Tl[) : 

dp(a) = 0.656 - (a ^° )3 '' — , a<1.15AU (41) 
W 1 + (a/a ) 3 - 6 a ^ ; 

where ao = 0.085 AU. We then generate an initial dis- 
tribution in the period ratio V (V > 1) of a two- 
planet system by generating two random variables, a\ 
and a>2, from this probability distribution and comput- 
ing V = max{(ai/a2) 3 / 2 , (a 2 /ai) 3 / 2 }. This procedure 
generates a smooth initial distribution in period ratio. 
For the numerical integrations, we set the semi-major 
axis of the massive planet to unity and the semi-major 
axis of the test particle is then V~ 2 ^ 3 . 

In some simulations the test particles have non-zero 
eccentricities e and/or inclinations i. These are assumed 
to be randomly distributed following a Rayleigh law, 

dp= ^exp(-ix 2 / ( 7 2 ), (42) 

where x = e or i and o~ x is an input parameter that is 
related to the mean and rms eccentricity or inclination by 
(x) = ^JV/2a x = 1.253cr x , (x 2 ) 1 / 2 = y/2a x = lAUa x . 

We treat the planet and host star as point particles, 
i.e., we do not account for possible collisions of the test 
particles with either body. 

5.2. Numerical results and comparison with analytic 

theory 

We start by considering two simple fiducial models to 
compare with our previous theoretical results. First, we 
consider a model for a slowly growing planet mass, which 
contains 10 4 test particles and has planet-formation 
timescale t p \ = 10 4 [27r/n p i] (or 10 4 orbits of the ex- 
terior planet). Second, we consider a model for a rapidly 
growing planet mass, which contains 2 x 10 4 test parti- 
cles and starts with the final planet mass, i.e., t p \ = 0. 
Both models have a final planet mass m p i = 0.001m* 
(1 Jupiter mass), a planet eccentricity e p i = 0, and zero 
eccentricities and inclinations for the test particles (i.e., 
o~e = o~i = in Equation I42|). The initial distribution 
of semi-major axes or periods of the test particles is ob- 
tained using the procedure in £|5.1[ with the period ratio 
restricted to the range 1.3 < V < 2.5. 

In the simulation for the slow case, ss 4% of the test 
particles were lost to escape orbits, mostly from orbits 
initially near the planet with Vi < 4/3, while for the fast 
case the same happens for rj 10% of the test particles, 
almost all from orbits with Vi < 7/5. 

5.2.1. Final eccentricity 

In Figure [5j we compare our numerical results for the 
final eccentricity as a function of the initial period ratio 
with those obtained from the analytic formalism of SQ] 
for a slowly growing mass planet. As expected, substan- 
tial eccentricities are excited in orbits that start close 
to the first-order (4:3, 3:2 and 2 : 1) resonances. The 
eccentricities are also excited in smaller ranges around 
the second-order (7:5 and 5:3) resonances, a result that 
is not captured by our simple analytic theory. Analytic 
curves for the 3 : 2 and 2 : 1 resonances (blue and red, 
respectively, from Equations [34] and [35]) are plotted for 



0.2r 

0.18 



0.16 




1.4 1.6 1.8 2 2.2 



initial period ratio 

Fig. 5. — Final eccentricity as a function of the initial period ra- 
tio ni/n p \ for a simulation with 10 4 test particles and parameters 
m pl = 0.001m,, t p i = 10 4 , e p i = 0, and <r e = 0. Results are shown 
after 5 X 10 4 planet orbits. The blue and red lines indicate the 
analytic results for a slowly growing planet, for the 3 : 2 and 2 : 1 
resonances, respectively. The analytic results are only shown for 
particles that do not cross the resonance separatrix, since parti- 
cles that do cross have librating eccentricities (these produce the 
elongated ovals filled with scattered points). The horizontal lines 
denote the approximate analytic prediction for the me an c cccntric- 
ity of particles with oscillating eccentricity (Equation 1330 . 

orbits that do not cross the resonant separatrix, since 
these have stationary forced eccentricity (compare the 
red dashed lines in Figure 0]). The analytic curves agree 
quite well with the simulations close to the resonances. 
The deviations grow as we move away from the reso- 
nances, since the approximation of a single resonance be- 
comes less and less accurate. For orbits that do cross the 
separatrix, the figure shows an approximate estimate for 
the mean eccentricity from Equation (j33| as a horizontal 
line; the extent of this line marks the range of period 
ratios corresponding to such orbits, from Equation (|36[) . 
In this snapshot of the simulation, the regions in which 
the forced eccentricity oscillates are marked by clouds of 
points; the widths of these clouds agree reasonably well 
with the analytic theory, although the theoretical esti- 
mates of the mean eccentricities, marked by horizontal 
red and blue bars, are somewhat too high. 

For the case when the planet mass grows quickly the fi- 
nal eccentricity is never stationary: it oscillates between 
zero and a maximum value similar to that in Figure [5] 
The results obtained from the resonance Hamiltonian in- 
tegration in SJ3] agree quite well with the numerical simu- 
lations. Near the 2 : 1 and 3 : 2 resonances, both produce 
a similar cloud of final eccentricities for a given initial 
period ratio. 

5.2.2. Period-ratio distribution 

In Figure [6] we show the distribution in period ratio for 
the test particles in the fiducial simulations. The upper 
and lower panels show the results for the slowly and fast 
growing mass cases, respectively. We also plot in red lines 
the probability density function (PDF) around the first- 
order resonances, as computed in $3]and £14.2l respcctivelv 
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TABLE 1 

Equivalent widths around resonances 



1.8 2 

period ratio 

Fig. 6. — Number density of test particles as a function of 
the period ratio V (black histograms), for the fiducial three- 
body integrations with parameters m p i = 0.001m*, e p i = 0, 
<r e = 0, and <t; = 0. The initial distribution in period ratio is 
given by Equation l|41jl and shown by blue lines. The upper 
panel shows the distribution after 5 X 10 4 planet years for a 
simulation with growth time t p \ = 10 4 planet years, as well as 
the theoretical PDF around the 4:3, 3:2, and 2 : 1 resonan ces, 
shown in red lines. The latter distribution is obtained as in i|4.2l 
(see Figure |4j and rescaled to the number of systems in the 
three-body integrations. The lower panel shows the same for a 
set of three-body integrations in which the growth time i p i = 0; 
in this panel the theoretical PDF is obtained as in f|3](sce Figure[3). 



(see black lines in figure [3] and 2]) , scaling our results to 
the mean number of initial planets in the simulation as 
expected from the background distribution in Equation 

As predicted by the analysis for a slowly growing planet 
in 21 there is a peak in the number of planetary systems 
at period ratios slightly larger than the exact resonances 
and a deficit at period ratios smaller than the resonance 
position. The strongest enhancements are at the first- 
order 4:3, 3:2, and 2:1 resonances, but the effect is 



Resonance 


EW+ 


EW_ 


|EW±| 
from Equation H39H 


4:3 


0.0073 a 


-0.034 


0.033 


7:5 


0.0043 


-0.0051 




3:2 


0.033 


-0.034 


0.033 


5:3 


0.0033 


-0.0036 




2:1 


0.030 


-0.030 


0.034 



Note. - Equivalent widths (Equation I29H from 
three-body integrations, for slow growth of the planet 
mass. 

a Note that 82% of the test particles around 4 : 3 are lost 
to escape orbits and this decreases the value of EW+, 
while slightly increasing the absolute value of EW_ . 

also weakly visible at the second-order 7 : 5 and 5 : 3 res- 
onances. Note, however, that the 7 : 5 resonance is close 
enough to 4 : 3 that the assumption of an isolated reso- 
nance used in the theory is suspect. Similarly, the 4:3 
resonance is close enough to the planet that a signifi- 
cant fraction of the particles at smaller period ratios are 
ejected by the planet. 

The features seen in the three-body integrations at 
the first-order resonances are roughly consistent with the 
predictions of 21 (red lines in Figure EJ)). Specifically, the 
theory is able to reproduce reasonably well the position 
and height of the peaks and the width and depth of the 
gaps at the 3 : 2 and 2 : 1 resonances. 

Similarly, as seen in the lower panel of Figure [6l the 
theory of fj3]for the rapidly growing planet agrees reason- 
ably well with the simulations around the 3 : 2 and 2 : 1 
resonances. We observe, however, that the gap narrow 
of the 3 : 2 resonance is deeper than the single-resonance 
theory predicts, which is mainly due to a fraction of the 
particles that are excited to high eccentricities (e ~ 0.2- 
0.6) and larger period ratios (310 particles with initial 
period ratio 1.4 < Vi < 1.5 have a final period ratio 
1.6 < V < 2.5, while 16 are ejected from the system). 
The agreement with the single-resonance approximation 
is even worse closer to the planet: in fact, this approxima- 
tion breaks down completely when first-order resonances 
overlap, which is ex pected for perio d ratios < 1.33 for a 
Jupiter-mass planet (jWisdomlll98Q[ ). Additional discrep- 
ancies arise from ejection of test particles by the planet, 
which occurs for 77% of the particles with Vi < 7/5. 

In summary, the single- resonance analysis of <J3] and 
is able to reproduce the main features of the three-body 
integrations around the 3 : 2 and 2 : 1 resonances. The 
period-ratio distributions for both slow and fast planet 
growth exhibit peaks wide of the resonance (period ra- 
tios larger than the resonance value) and troughs or gaps 
narrow of the resonance. This result suggests that a 
broad range of planet-formation processes on intermedi- 
ate timescales will produce a similar period-ratio distri- 
bution, so long as the planet does not migrate. However, 
the distributions differ in detail when the planet mass is 
increased slowly or fast. Specifically, the slow case pro- 
duces both peaks and troughs that are sharper and lie 
closer to the resonance. 

5.2.3. Equivalent widths around resonances 

The equivalent widths measured from the fiducial simu- 
lation of a slowly growing planet are shown in Table [1] 
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along with the predictions from Equation (|39l) . In most 
cases the EWs on either side of the resonances are almost 
equal and opposite, confirming that the resonances lead 
to shuffling of the test particle semi-major axes but not 
any overall loss or accumulation. 

As expected, the largest EWs are at the 3 : 2 and 2 : 1 
resonances, ~ 0.033 and ~ 0.030, respectively. The EWs 
at 3 : 2 and 2 : 1 agree with the single-resonance approx- 
imations (Equation [29]) . The 4:3 resonance agrees less 
well because 182 test particles around this resonance were 
ejected. But even in this case the single-resonance theory 
matches at least the width of the trough, EW_. 

For the fast mass increase case, we measure EW_ = 
-0.0262 and EW+ = 0.0254 around the 2 : 1 resonance, 
and the single-resonance theory (Equation I30|) yields 
|EW±| = 0.0244, in reasonably good agreement. At the 
3 : 2 resonance, EW_ = -0.037 and EW+ = 0.030, while 
the single-resonance theory predicts |EW±| = 0.0237. 
Here the difference between the single-resonance theory 
and the three-body integrations is mostly due to the ex- 
citation of somewhat larger eccentricities in the integra- 
tions that tend to populate the peak wide of 3 : 2 with 
more particles, while excavating a deeper gap around this 
resonance; 310 test particles with 1.4 < V% < 1.5 are pro- 
moted to orbits with V > 1.6 and 16 to escape orbits, 
which accounts for the difference in magnitude between 
EW_ and EW+. We have checked that reducing the 
mass of the perturber improves the agreement between 
the single-resonance theory and the integrations. 

5.3. Dependence on input parameters 

All of the simulations described so far consider a Jupiter- 
mass planet with zero eccentricity and test particles 
in initially circular and coplanar orbits, so that a di- 
rect comparison with the single-resonance theory can be 
made. Here we comment on how the distributions in pe- 
riod ratio are affected by different values of the planet 
mass and eccentricity, and the initial rms eccentricity 
and inclination of the test particles. 

5.3.1. Varying the planet mass m p \ 

We have run sets of three-body integrations in which the 
mass of the planet varies over the range 10 -4 -10~ 3 m». 
These experiments show that the integration results obey 
the scalings predicted by the single-resonance theory: 
in particular the eccentricity of the test particles scales 

as s e cx ™p{ 3 , while the equivalent width varies as 

|EW±| cx s cx Wpf 3 . Thus, varying the planet mass sim- 
ply produces re-scaled versions of the PDF around the 
3 : 2 and 2 : 1 resonances shown in Figure [6] 

These scalings are more difficult to study for reso- 
nances that are closer to the planet, both because of 
particle ejections and possible resonance overlap. Nev- 
ertheless, in most cases the analytical single-resonance 
theory allows us to scale predictions for the distribution 
of period ratios to different resonances and different plan- 
etary masses. 

5.3.2. Varying the eccentricities: e p i and a e 

Varying the eccentricities of the test particles can 
strongly modify the final distribution of period ratios. 
In particular, values of er e (Equation [42]) that arc larger 



than the eccentricity scale s e (Equation ll9[) tend to sup- 
press the characteristic peak+trough or "P-Cygni" fea- 
ture around first-order resonances. For reference, from 
Equation ([T9]) the eccentricity scale for a Jupiter-mass 
planet at the 3 : 2 and 2 : 1 resonances is 0.05 and 0.063, 
respectively. For values of a e lower than s e the PDFs 
are similar to those starting with initially circular orbits 
(a e = 0). 

Non-zero eccentricity of the perturber is expected to 
introduce chaotic zones in the vicinity of the resonances, 
such that initially circular orbits are excited to higher ec- 
centricity over a wider range in period ratios compared 
with the case of a zero-eccentricity perturber. In numer- 
ical simulations with perturber eccentricities e p i = 0.05- 
0.1, we observe that a large fraction of test particles with 
Vi < 3/2 are ejected to escape orbits, leaving a strong 
gap around 3 : 2 and almost no evidence of a peak wide 
of the resonance. Additionally, the gap around the 2:1 
resonance is observed to be wider and less deep than 
for a zero-eccentricity perturber, while the peak wide of 
the resonance is observed to be somewhat weaker. Non- 
zero eccentricities could also lead to slow evolution of the 
period-ratio distribution on timcscales much longer than 
the 10 4 planet orbits used in the fiducial simulations; 
we have carried out extended integrations for up to 10 5 
planet orbits but see no additional changes in the period- 
ratio distribution. In summary, the non-zero eccentric- 
ity of the perturber tends to smear out the characteristic 
"P-Cygni" profile that we observe in the zero-eccentricity 
simulations. 

5.3.3. Varying the inclinations: o~i 

We have run three-body integrations with parameters 
<j e = 0.01 and <x; = 0.1. The latter parameter gives a 
mean inclination of (i) w 7°, a value that is somewhat 
higher than estimates of the mean inclinations in Kep ler 
multi-pla net systems: < 5° (|Tremaine k, Dondl2011|) or 
1.0°-2.3° (jFabrvckv et al.ll2012[) . Although the dominant 
features around the first-order 3 : 2 and 2 : 1 resonances 
are still present, the "P-Cygni" profile shows a less sharp 
peak that is shifted to larger period ratios, and gaps that 
are slightly wider than in the zero-inclination case. More 
specifically, we find that at 2 : 1 the equivalent widths 
slightly increase to |EW±| = 0.033 (see values in Table 
[T]for comparison with the fiducial simulation), while at 
3 : 2 we observe that more particles are ejected to escape 
orbits relative to the zero- inclination case, resulting in an 
increase in magnitude of EW_ to —0.038 and a decrease 
of EW + to 0.021. Finally, the gap at the second-order 
7 : 5 resonance gets wider by a factor of roughly two. 

6. COMPARISON WITH THE KEPLER SYSTEMS 

We now compare our results with the Kepler catalog of 
242 tw o-planet systems as obtained from lFabrvckv et al.1 
(]2012|) . We do not consider systems with more than two 
planets since we have not investigated three-body res- 
onances in this paper. We also compare with the 33 
two-planet systems discovered by radial velocity (RV) 
observations, using t he The Exoplanet O rbit Database as 
of September 2012 ([Wright et all 120111 ). We limit our 
analysis to a range in period ratio V of 1.3-2.5, in which 
the first-order resonances are 2:1, 3:2, and 4:3. After 
this cut our sample is reduced to 116 and 10 two-planet 
systems from the Kepler and RV catalogs, respectively. 
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period ratio period ratio period ratio 

Fig. 7. — Number density of systems as a function of the period ratio. The left panel shows the Kepler sample (black histogram) and 
the Kepler sample plus the RV sample (blue line) . The middle and right panels show the results of the fiducial simulations for a slowly 
and rapidly growing planet, respectively. The solid red lines indicate the initial distribution in period ratio given by Equation H41| l. scaled 
to the number of systems in each panel (ignoring the small RV contribution in the left panel) and the red dashed lines its 1— a error bands. 
The vertical dot-dashed lines show the positions of the first- and second-order resonances. All histograms have the same binning of 0.05. 




-0.04 -0.02 0.02 0.04 -0.04 -0.02 0.02 0.04 -0.04 -0.02 0.02 0.04 



V/Vres ~ 1 V/V res - 1 P/Vre. - 1 

Fig. 8. — Distribution in period ratio of the Kepler two-planet systems near the 3:2 and 2: 1 resonances (left and middle panels). In 
the right panel the two distributions are added together (for given values of the shifted and normalized period ratio V fPres — 1)- The red 
line indicates the theoretical distribution for a slowly growing planet (the same as in Figure [4{ with a strength parameter (Equation I18H 
s = 0.01, normalized to the number of Kepler planets in each panel. 



In Figure wc show the distribution of period ratios 
for the observations (left panel) and the fiducial simula- 
tions. We also show the expected initial distribution of 
period ratios as determined by the algorithm described 
in the paragraph containing Equation (|41[) , along with 
its 1-er confidence limits to suggest whether peaks and 
dips are significant. 

The observations exhibit statistically significant peaks 
just outside (larger V than) the 3 : 2 and 2 : 1 resonances, 
although the peak at 2 : 1 only exceeds 1.5-cr significance 
when the RV planets are included in the sample (these 
statements depend on the binning, which has been cho- 
sen to maximize the significance of the near- resonant fea- 



tures). The observations show troughs narrow of the 3 : 2 
and 2 : 1 resonance, although these are only marginally 
significant (1-1. 5a). There is also a deficit of systems 
with period ratios V < 1.4; this likely arises from the 
depletion of planets that suffer close encounters and are 
scattered onto collision or ejection orbits. There are no 
statistically significant features at any second-order res- 
onance. Therefore, we will concentrate our analysis on 
the first-order 2 : 1 and 3 : 2 resonances. 

6.1. The 3:2 and 2:1 resonances 

As pointed out by iLissauer et al.l ()2011b[ ) and 
iFabrvckv etldl $201% and also shown in Figure [3 
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two-planet systems close to the 3 : 2 and 2 : 1 resonances 
prefer period ratios wide of the resonance within a few 
percent; in particular, the data show significant peaks 
wide of the resonances and perhaps a trough narrow of 
the resonance in the case of the 2 : 1 resonance. Here we 
compare the distribution of planetary systems around 
the 3 : 2 and 2 : 1 resonances with our simple models to 
test whether they are consistent. 

Hereafter we shall concentrate on the more realistic 
model in which the mass grows slowly and will refer to 
the distribution obtained using the single-resonance ap- 
proximation (2] and Figure [4|) as our theoretical PDF. 

In Figure [5] we show the distribution of systems as 
a function of the shifted and normalized period ratio 
V /Tics — 1. Positive and negative values of this variable 
lie wide or narrow of the resonance, respectively. For 
comparison, we show our theoretical PDF normalized to 
the number of planets in each panel. Here we just con- 
sider one value of the strength, s = 0.01, which is equiv- 
alent to a planet mass m p i ~ 0.3A'/,/ or m p \ ~ OAMj at 
3 : 2 and 2:1, respectively. The corresponding equivalent 
width is EW± = ±0.00956 7> rcs (see Equation [3$]). 

We observe in the left panel of Figure [8] that the the- 
oretical PDF produces a peak similar in height and lo- 
cation to that observed in the Kepler sample around the 
3 : 2 resonance, but the Kepler sample shows no sign of 
the expected gap. In contrast, in the middle panel the 
gap inside the 2 : 1 resonance seems to be larger in the 
Kepler observations than in the theoretical PDF. Finally, 
in the right panel where we add the 3 : 2 and 2 : 1 reso- 
nances together, the overall shape of the period distri- 
bution seems to agree, at least visually, with that of the 
theoretical PDF. More quantitative comparison requires 
statistical tests, to which we now turn. 

6.1.1. K-S test 

We use a Kolmogorov-Smirnov test (K-S test) to deter- 
mine the probability that the observed sample is drawn 
from the same distribution as our theoretical distribu- 
tion. Since we have found that the PDF obtained from 
the single-resonance model is a good approximation to 
the PDF obtained from three-body integrations starting 
with low-eccentricity, low-inclination orbits, we shall use 
the single-resonance predictions for the subsequent anal- 
ysis. 

Note that the Kepler sample contains planets with a 
wide range of masses. Thus we are fitting the data to 
some loosely defined "typical" value of the planet mass 
mpi or dimensionlcss resonance strength s (Equation llSp . 

In Figure we show the p-values obtained from the 
K-S test for the Kepler planets near the 3 : 2 and 2 : 1 
resonances, and for the sum of these two distributions at 
given V/Vrcs — 1 (from top to bottom panels). For ref- 
erence, we also show the p-value obtained by comparing 
the Kepler planets to a uniform distribution (horizontal 
dashed line). In this test, we have restricted the samples 
to the range \V /V rC s — 1| < 0.04. We have experimented 
with other ranges between 0.03 and 0.1 and found that 
0.04 provides the most stringent comparisons. For ranges 
larger than 0.04 there are too many particles far from 
resonance that dilute the signal; also the test becomes 
biased by large-scale gradients in the period ratio distri- 
bution. For ranges smaller than 0.03 there are too few 
planets left in the sample. 




Fig. 9. — p-values as a function of the resonance strength s 
(Equation 1181 obtained from a Kolmogorov-Smirnov test compar- 
ing the Kepler two-planet sample with the theoretical PDF for a 
slowly growing planet that does not migrate (solid line). The hor- 
izontal dashed line represents a uniform density distribution, and 
the dash-dot line represents a significance level p = 0.05 — models 
below this line are excluded at the 95% confidence level. We limit 
the test to the range IP/Prcs — 1| < 0.04 to maximize its power. 

As usual, if the p— value is below a given value a there 
is a probability of 1 — a that the Kepler data are not 
drawn from the PDF predicted by the single-resonance 
model. Thus, in each panel we add a horizontal dot- 
dashed line at 0.05 to discriminate between models at 
the 95% confidence level. The middle panel of Figure [8] 
shows that there is no statistically significant signal at 
the 2 : 1 resonance, in the sense that a uniform distri- 
bution of period ratios is consistent with the data (at 
p = 0.53). This result does not exclude the possibility 
that a significant signal would be revealed by some other 
more sensitive test. In the top and bottom panels, there 
is a significant signal in that the uniform distribution is 
excluded at about the 97% confidence leyel. 

The se results differ from those of iFabrvckv et al.l 
(|2012f ). who find significant evidence that the distribu- 
tion in period ratios around the 3 : 2 and 2 : 1 resonances 
are not drawn from a smooth distribution (p = 0.0046 
and 0.029, respectively), both more than an order of mag- 
nitude smaller that our p-values. We attribute this dif- 
ference mainly to our restricted data set: we use only 
two-planet systems whereas Fabrycky et al. use all planet 
pairs in multi-planet systems, so our sample is ~ 3 
times smaller. The difference might also hint that near- 
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resonant features are more common in multi-planet sys- 
tems than in two-planet systems. 

7. DISCUSSION 

The few hundred multiple planet candidates discovered 
by the Kepler spacecraft provide a unique window into 
how planets form and evolve, although we do not yet 
know how to interpret the limited view that we have 
through this window. This paper is focused on the distri- 
bution of planets near mean-motion resonances in multi- 
planet systems, which has two important features. 

First, the number of systems near mean- motion reso- 
nances is small: the left panel of Figure [7] shows that the 
only resonance feature that is significant at more than 
the 1.5-cr level is the peak at the 3 : 2 resonance. As dis- 
cussed in SJTJ this result is striking because since planet 
migration is believed to be a common process, and cap- 
ture into a first-order resonance during convergent mi- 
gration is certain if the eccentricities and inclinations of 
the planet pair are sufficiently small and migration is 
sufficiently slow when the resonance is crossed. 

Second, at both the 3 : 2 and 2 : 1 first-order resonances 
the peak in the period-ratio distribution appears just 
wide of the resonance (i.e., at period ratios greater than 
the resonance value, if the period ratio is defined to be 
> 1), rather than at the resonance. In addition, the 2 : 1 
resonance appears to have a trough in the period-ratio 
distribution just narrow of the resonance, although the 
statistical significance of the trough is w eak. All of these 
featur es were pointed out already by iFabrvckv et al.l 

dmJ). 

Given these findings, it is natural to ask what distribu- 
tion of period ratios should be expected near resonance 
if there is no migration, i.e., if planets form in situ. We 
have investigated two extreme limits of this process, in 
which the planet forms fast or slowly (relative to the 
characteristic dynamical time associated with the reso- 
nance) . We find that in both cases the distribution of test 
particles near resonance develops a shape characterized 
by a peak wide of the resonance (period V > V res where 
the period ratio at resonance "P res > 1) and a trough 
narrow of the resonance (see Figures [3] and 2} . The peak 
and trough result from a redistribution of particles from 
period ratios V < V lcs to V > 'Pros, rather than from an 
overall gain or loss of particles as would occur through 
processes such as ejection or migration. 

7.1. Planet masses in the Kepler sample 

The mean and media n radius for the planet s in the 242 
two-planet systems in IFabrvckv et all (|2012| ) are 2.46i?© 
and 2.17i?0, respectively To compare our models to 
the data, i.e., to estimate the strength parameter s in 
Figure IH1 we need the mass-radius relation for the Kepler 
planets. 

(i) Fitting to the planets in the solar system, 
IFabrvckv et all (12011) find M = M e (i?/i? e ) 2 06 for 
R > i? e , which yields M = 6.4M e and 4.9M ffi for the 
mean and median mass in the Kepler two-planet sam- 
ple, (ii) From a statistical fit to transit timing v aria- 
tions in the Kepler sample, IWu fc Lithwickl (|2012l ) find 
M = 3M e (i?/i? e ), which implies 7.4M e and 6.5M e for 
the mean and median mass, (iii) There are four planets 
with measurements of both mass and radius in the ra- 
dius range 2-3 i?©; these have mean and median masses of 



7.0Mq and 6.1Mq. All three of these crude approaches 
suggest a typical planet mass of about 6-7Mq for the 
Kepler planets. 

7.2. Planet masses required to explain features in the 
period distribution near resonances 

Our simple theoretical model, based on slow growth of 
the planet masses with no migration, is consistent with 
the data (p > 0.1) shown in Figure [9] for values of the 
strength s between 0.006 and 0.0180. Equation (fT8|) gives 
s = 0.023(m pl /A/j) 2 / 3 and s = 0.0179(m pl /Afj) 2 / 3 at 
3 : 2 and 2:1, respectively; focusing on the 3 : 2 resonance, 
which contains the strongest signal, consistency then re- 
quires 40 < m p i/M® < 220. Our calculations are for a 
system containing a test particle and a planet of mass 
m p i; in the Kepler two-planet systems both planets have 
non-zero mass and in the absence of a detailed numerical 
study of this case the most appropriate choice is proba- 
bly to identify m p i in these formulae with twice the mean 
planet mass. Thus the mean planet mass required to ex- 
plain the resonant structure seen in the Kepler data at 
the 3 : 2 resonance is roughly 2O-1OOM0. 

The low end of this mass range is about three times the 
mean planet mass estimated in the preceding subsection. 
There are several possibilities for bridging this gap: 

Test-particle calculations underestimate the dynamical effect 
of resonances — Our analytic calculations are only pos- 
sible because one of the two planets was approximated 
as a test particle. If both planets have non-zero masses, 
the dynamical problem acquires extra degrees of free- 
dom; even the planar single-resonance model has two 
degrees of freedom which allows chaos caused by over- 
lap of the closely spaced resonances with critical angles 
(p + l)Aa — pXi — LO\ and (p + 1)A2 — pAi — u>2- We have 
carried out exploratory integrations in which we replace 
the three-body system containing a planet of mass m pl 
and a test particle with two planets of mass m p \/2, leav- 
ing the other initial conditions and parameters the same. 
We find that the peaks and troughs at the 3 : 2 and 
2 : 1 resonances become significantly stronger, and that 
a secondary peak develops narrow of the 2 : 1 resonance. 

Long-term evolution — Systems with more than two de- 
grees of freedom can exhibit slow evolution due to 
weak chaos (a classic example is the Kirkwood gaps in 
the asteroid belt , which can evolve on Gyr timescales; 
iMorbideTTil 119961 ) . Integrations that follow two plan- 
ets with non-zero masses, eccentricities, and inclinations 
may have resonant features that evolve over timescales 
much longer than those we have examined here. 

Underestimated planet masses — Simultaneous measure- 
ments of mass and radius are available for 8 planets in the 
radius range l-3i?0. These show an order of magnitude 
range of mean density, from 0.7^°'^ g cm~ 3 for Kepler 
llf (jLissauer et all 1201 laf ) to 8.8+2'g g cm -3 for Kepler 
10b (jBatalha et al.ll2011h . This wide ran ge implies either 
a diverse set of planetary compositions (|Wu fc Lithwickl 

6 We recognize that the KS test is designed for hypothesis test- 
ing, not parameter fitting, but the data and the models do not 
justify a more sophisticated approach. 
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l2012f ) or unmodeled errors in some or all of the measure- 
ments. In either case the actual value of the "typical" 
Kepler planet mass is far more uncertain than the sim- 
ple estimate of 6-7M® obtained above. In particular, if 
there is a substantial population of rock-iron planets with 
radii of 2-3i?0 these would have much la rger masses: for 
a radius of 2.3i?®, iSeager et al.1 (|2007t ) find masses of 
30-200M m f or a v ariety of rock-iron compositions, and 
ISwift et al.1 (|2012h find' masses from 30M§ for basalt to 
lOOAf^ for nickel-iron planets. 

7.3. Relation to previous work 

IBatvgin fc Morbidellil (pfillh and lLithwick fc Wul (pOllh 
have shown that tidal dissipation can repel planet pairs 
from exact resonance, producing a peak wide of the reso- 
nance and a trough narrow of the resonance as observed 
in the Kepler data. However, this process requires that 
(i) the planet pairs form in resonance; (ii) the correct 
amount of dissipation is present to repel the planet pairs 
from resonance by a percent or so. 

Moreover, if tidal dissipation from the host star were 
responsible for shaping the period-ratio distribution 
around resonances, one should expect a strong depen- 
dence on orbital radius or period (the characteristic 
timescale for changes in semi -major axi s a varies as a 8 in 
the ti dal evolution model of lHutllT981l ). However. iReinl 
(|2012D points out that the period-ratio distribution in the 
Kepler sample restricted to inner-planet periods < 5 d 
looks identical to the distribution restricted to inner- 
planet periods > 5d. We have repeated Rein's analy- 
sis for the restricted sample of two-planet systems used 
here and observe that there is a hint that the resonant 
features are more pronounced at the 3 : 2 resonance for 
the longer period planets, and more pronounced at the 
2 : 1 resonance for the shorter periods planets. However, 
the number of planet pairs in our restricted sample is too 
small to make a statistically significant detection of any 
d ifference. 

IReinl (|2012f ) has also investigated whether stochastic 
migration can reproduce the Kepler period-ratio distri- 
bution. He argues that conventional "smooth" migration 
cannot reproduce the observations because it produces 
too many planet systems in exact resonance. However, 
these large pile-ups at resonance are smeared out by in- 
cluding stochastic forces that might be expected from 
the likely turbulent nature of the protoplanctary disk. 
He shows that for the correct combination of stochas- 
tic and smooth migration forces, the final period-ratio 
distribution looks similar to that of the Kepler planets. 

We note that eccentricity measurements of near- 
resonant planets do not distinguish between these mech- 
anisms: tidal dissipation damps the free eccentricity so 
the residual eccentricity is equal to the forced eccentric- 
ity, which is just the fixed point of the single-resonance 
Hamiltonian shown in Figure [5] Slow growth of the 



planet causes a test particle on an initially circular orbit 
to follow the fixed point as the planet mass grows, so 
that it will also have zero free eccentricity unless it has 
crossed the separatrix. 

8. SUMMARY 

We have studied the orbital distribution of two-planet 
systems near first-order mean-motion resonances in the 
simplest possible model of planet formation: there is no 
energy dissipation or migration, and planets form in situ 
on circular, coplanar orbits, with masses growing at a 
prescribed rate. We have examined whether this toy 
model can explain the signatures in the period-ratio dis- 
tribution near resonances that are observed in the Kepler 
sample of 242 two-planet systems. 

Our approach is to construct a simplified Hamiltonian 
that isolates the perturbations due to a given first-order 
resonance. We then solve for the long-term dynamics of a 
test particle in this Hamiltonian, focusing on two limiting 
cases: rapid and slow mass growth of the planet. We 
have used numerical integrations of the restricted three- 
body problem to confirm that the resonance Hamiltonian 
captures the main features of the orbit evolution. 

We find that the distribution in period ratios resembles 
a "P-Cygni" profile, where orbits are evacuated from nar- 
row of the resonance (i.e., closer to the perturbing planet) 
and pile up in regions wide of the resonance. These fea- 
tures are present whether the planet grows fast or slowly, 
though they are stronger in the latter case, and hence 
should be present for a wide range of growth histories. 
The resulting structure in the period-ratio distribution is 
strongly reminiscent of the peak-trough feature observed 
at the 3 : 2 and 2 : 1 resonances in the Kepler systems. 

The model and observed period-ratio distributions 
near these resonances are consistent for mean planet 
masses in the range 2O-1OOM0. This is larger than ex- 
pected from the handful of Kepler planets with typical 
radii of 2-3i? ffl and measured masses by a factor of 3-15. 
However, the mass-radius relation for Kepler planets is 
quite uncertain and shows large scatter, and the required 
masses are consistent with the masses expected for plan- 
ets in this radius range that are composed of a mixture 
of silicates and iron with no extended atmosphere. Othcr 
effects such as weak chaos operating on Gyr timescalcs 
may also narrow the gap between the model and the 
observations. Our model suggests that the resonant fea- 
tures seen in the Kepler multi-planet systems may not 
require either dissipation or migration during the planet- 
formation process. 

We are grateful to Subo Dong, Darin Ragozzinc, and 
Hanno Rein for their insights and suggestions. RM grate- 
fully acknowledges support from the Institute for Ad- 
vanced Study, and CP acknowledges support from the 
CONICYT Biccentenial Becas Chile fellowship. 



APPENDIX 

MIGRATION SPEED AND CAPTURE INTO RESONANCE 

Capture into a (p + 1) : p resonance during convergent migration is certain if the planets cross the resonance slowly 
enough and their eccentricities and inclinations are small enough. These statements can be quantified using the 
formalism of Sj2j 

In convergent migration, the dimcnsionless resonance distance A (eq. 1 1 Tj> is increasing with time. Since the resonant 
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Hamiltonian K (eq. 1161) is dimcnsionless, resonance crossing without capture should occur when dA/dr > g where g 
is some constant of order unity. Converting this inequality to physical units yields 



d 



logn p i>^|p 2 (p+l)| 1 /8„ pl 



TOpl 



fp 



4/3 



(Al) 



Define the migration time of the massive planet to be t m - lg = (d\ogn p \/dt) , then the condition for resonance crossing 
without capture is 

3.0 x 10 4 yr P pl /10M ffi m* x 4/3 



100 d V m p i M Q 



(A2) 



where P p \ is the massive planet's orbital period and 



\p 2 {p- 



1 resonance) and 
ic 



T - l)|^^ ^p^ is 1.36 3 for — 1 (2 . 
3.377 for p = 2 (3 : 2 r esonance ) . The constant g is estimated to be 2.5 bv lFriedfa nd (200l|) from approximate analyt 
calculations and 2.7 bv iQuillenl ()2006l ) from numerical orbit integrations. 

The maximum eccent ricity at which capture into resonance is certain during convergent migration has been derived 
by several authors (e.g.. lHenrard fc Lema itrc 1 ?>'-•>: iBorderies fc Gold reich 19841). If the eccentricity of the test particle 
when far from resonance is eo, then the area enclosed by its orbit in the x-y plane is Aq = ite\j s\ (eq. IT9]) . Because 
of adiabatic invariance, this area is conserved as the planet migrates. The test-particle orbit will be interior to the 
homoclinic orbit when it first appears, at A = 1 (cf. Figure [I), if A) < 67r, the area of the homoclinic orbit. As the 
planet continues to migrate and A continues to grow, the area inside the homoclinic orbit grows so the test particle 
remains inside the homoclinic orbit, i.e., it remains captured in the resonance. Thus the condition that resonance 
capture is certain is eo < V6s e or 



eo < 



2 i/2 3 i/e 



TOpl , 

Jp 

TO* 



1/3 



|p2(p+l)|Vfl 

where h p = \p 2 (l + p)\~ 2/9 \f P \ 1/3 is 0.909 for p = 1 and 0.728 for p = 2 
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